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Abstract. We study heat transport in a class of stochastic energy exchange systems 
that characterize the interactions of networks of locally trapped hard spheres under 
the assumption that neighbouring particles undergo rare binary collisions. Our results 
provide an extension to three-dimensional dynamics of previous ones applying to the 
dynamics of confined two-dimensional hard disks [Gaspard P & Gilbert T On the 
derivation of Fourier's law in stochastic energy exchange systems J Stat Mech (2008) 
PI 1021]. It is remarkable that the heat conductivity is here again given by the 
frequency of energy exchanges. Moreover the expression of the stochastic kernel which 
specifies the energy exchange dynamics is simpler in this case and therefore allows for 
faster and more extensive numerical computations. 



Submitted to: J. Stat. Mech. 
1. Introduction 

The derivation of Fourier's law of heat transfer in insulating solid materials is a difficult 
problem which has been challenging theoretical physicists for close to two centuries pQ . 
However, recent works [2j [3] on models of confined particles in interaction have shed new 
light on this problem, setting the stage for a systematic derivation of Fourier's law from 
first principles. These works suggest that there is hope to achieve such a first principles 
characterization of heat transfer and prove the validity of Fourier's law in a particular 
class of insulating materials known as aerogels. 

In order to achieve such a characterization, the authors of [21 [3] have focused their 
attention on classes of models which, like aerogels, combine the collisional dynamics 
of gases with the spatial structure of solids. Starting from a Hamiltonian description, 
it was shown that the heat conductivity of such models is universally given by the 
frequency of collisions between gas particles. This universality manifests itself when the 
gas particles are individually trapped in a porous solid material and only rarely collide 
with each other, mostly rattling around their traps. Under this assumption, individual 
particles typically achieve local equilibrium states at their respective kinetic energies, 
ergodically exploring their trapping cells, before energy exchanges proceed. 
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This local equilibration mechanism is rather similar to the trapping mechanism of 
tracer particles in a periodic Lorentz gas which lead Machta and Zwanzig jl] to infer a 
stochastic approximation of mass transport in the Lorentz gas. In this approximation, 
the diffusion coefficient is identified, up to dimensional factors, with the rate of jump of 
tracer particles from cell to cell. 

Likewise, in our systems, this local equilibration naturally yields a stochastic 
description of the time evolution of the probability distribution of the local energies 
in terms of a master equation. Such derivation was extensively studied in |6], 
relating to the dynamics of confined two-dimensional hard-disks. Irrespective of the 
dimensionality of the underlying dynamics, the main difficulty when analyzing the 
transport properties of such stochastic systems is that, unlike with the Machta-Zwanzig 
models which deal with independent tracer dynamics, the energy distributions cannot 
be reduced to single cell distributions. Nevertheless, the transport coefficient -in this 
case the heat conductivity- is identified, up to dimensional factors, as the rate of energy 
exchanges. 

The purpose of this paper is to provide an extension of the results presented in 
[B] to the dynamics of trapped three-dimensional hard spheres. Starting from the 
Liouville equation which describes the time evolution of phase-space distributions of 
such systems, we consider the reduction of this equation to a stochastic evolution for 
the local energies under the assumption that collisions between neighbouring particles 
are rare compared to wall collisions. We subsequently show that the stochastic kernel 
which characterizes the energy exchanges between neighbouring cells has the symmetries 
described in [Bj which ensure the identity between the heat conductivity and frequency of 
energy exchanges. Furthermore, the expression of the stochastic kernel is much simpler 
in the case of three-dimensional hard spheres than in that of two-dimensional hard disks. 
The model is thus amenable to higher precision numerical simulations, which allows us 
to confirm the preceding arguments to a higher precision than had been previously 
obtained in the framework of underlying two-dimensional dynamics. 

The paper is organized as follows. A summary of the results described in [B] 
is presented in section [2j In section [3] we introduce mechanical models of periodic 
networks of confined hard spheres and discuss the necessary assumptions upon which 
these systems are amenable to a stochastic reduction. The statistical evolution of these 
systems is considered in section [4] and its reduction to a stochastic equation in section 
[5j with a detailed derivation of the stochastic kernel. The identity between the energy 
exchange frequency and the thermal conductivity is established in section |6j Numerical 
results supporting our theoretical arguments are presented in section [7j Conclusions 
and perspectives are drawn in section |HJ 

2. Summary of the main results 

Consider a system of N energy cells ex, ... ,6^, with stochastic exchange of energy among 
pairs of neighbouring cells. We assume that the statistical evolution is described by the 
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following master equation, 
<9t-P/v(ei, . . • , Cat, t) 

" N f r 

5^ / d?7[W / (e a + ?7,e;,-r7|e a ,e fe )PA r (. . . , e a + 77, ...,65-77, •••>*) 



N 

2 



a,b=l ' 



-VK(e a , e 6 |e a - r),e b + r])P N (. . . ,e a , . . . ,e b , ...,t) , (1) 

where W^(e a , 6b|e a — 77, e& + 77) is the stochastic kernel specifying the process of exchange 
of energy 77 between two neighbouring cells a and b at respective energies e a and e b , with 
— £6 < ?7 < e a - We will be concerned here with kernels that do not depend on the specific 
pair a, b of neighbouring cells. 

Given the two energies e a and e b , the characteristic time scale of energy exchanges 
between neighbouring cells a and b is determined by the collision frequency 

v(e a , e b ) = J di]W(e a , e b \e a - 77, e b + 77) , (2) 

whose canonical equilibrium average at temperature T we denote vb(T) = (v(e a , e b ))r ~ 

Vt. 

In the proper hydrodynamic scaling limit, this master equation yields the time 
evolution of the local temperatures, 

2 f N 

T a (t) = - I n de { e a P N (ei, ...,e N ,t), (3) 
a J i=1 

where d is the dimension of the underlying dynamics. This evolution turns out to be 
given according to Fourier's law, namely 

d t T{x,t) = -d x [K{T)d x T{x,t)\, (4) 

where the heat conductivity k(T) is 

k(T) = MT)- (5) 

The result ^ establishes an identity between a macroscopic quantity, the heat 
conductivity, and a microscopic one, the frequency of energy exchanges between two 
neighbouring cells. As shown in [B] , the derivation of this identity follows from equation 
([I]) in two independent ways, which rely on special symmetries of the kernel. 

These symmetries concern the equilibrium averages of the first and second moments 
of the energy exchanges, 

j(e a ,e b ) = J dr]r]W(e a ,e b \e a - r],e b + 77) , (6) 
h(e a , e b ) = J dr] r] 2 W(e a , e b \e a - 77, e b + 77) . (7) 

In [6], we showed that the kernel W associated with the confining dynamics of two- 
dimensional hard disks has the symmetries: 

(v{e a , e&)) T = ^ (i e a - e b )j(e a , e b )) T = ]- (h(e a , e b )) T . (8) 
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Given the master equation pi), there are two ways of computing the heat 
conductivity. The first one is to assume a non-equilibrium stationary state resulting 
from a temperature gradient. One then finds that the heat conductivity has expression 

K i T ) = g (( e <> _ e b)j(e a , e b )) T , (9) 

which, using the properties of the kernel ([8]), in its turn yields (|5]). 

The second way of computing the heat conductivity is through the Green-Kubo 
formula, from which it turns out that only static correlations contribute to the transport 
coefficient with expression 

K(T) = ±(h(e a ,e b )) T . (10) 

And, again using the symmetries of the kernel ([§]), we obtain (|5]). 

The goal of this paper is to provide a derivation of the master equation ([I]) associated 
with the stochastic energy exchanges of rarely interacting confined hard spheres and 
show that the corresponding kernel has the symmetries (pi). The identity (JsJ) ensues. 

We show in this paper that these results extend verbatim to the kernel associated 
with the confining dynamics of three-dimensional hard spheres. 



3. Networks of confined hard spheres 

Consider a lattice of confining three-dimensional cells, each containing a single hard- 
sphere particle. The mechanism of confinement prevents mass transport. However 
we imagine a form of semi-porosity by which particles in neighbouring cells are able to 
perform elastic collisions with each other, thereby exchanging energy. Such a mechanical 
system with hard-core confinement is depicted in figure [TJ Two dimensional versions of 
similar systems were extensively studied in [2J, [5] . We note that the geometrical details of 
the confining cells are irrelevant for what follows so long as the local dynamics mixes the 
velocity angles. That is, each cell taken individually with a single moving particle inside 
them is a semi-dispersing billiard with uniform equilibrium measure on the constant 
kinetic energy surface. This assumption validates the local equilibrium distributions 
described in the next section. 

These systems combine the collisional dynamics of gases and the spatial structure 
of solids. There is therefore a natural distinction between the local dynamics, which deal 
with the interactions between the moving particles and the solid matrix of their confining 
cell^fj and the interacting dynamics, by which two moving particles in neighbouring cells 
perform an elastic collision. On the one hand, the local dynamics are characterized by 
a wall-collision frequency, i^w, which depends on the geometry of the confining cell as 
well as on the kinetic energy of the moving particles through well known results of 
ergodic theory [7j. On the other hand, the interacting dynamics are characterized by 
the frequency of binary collisions, i.e. the frequency of collisions between neighbouring 
particles, denoted z/ B . 

| In our models, no energy exchanges take place between the confining walls and the moving particles. 
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Figure 1. Hard-sphere particle trapped in a cuboid cell with cylindrical edges. One 
imagines a material consisting of many copies of such cells forming a spatially periodic 
structure. The cells are semi-porous in the sense that particles are prevented from 
escaping, and can yet partially penetrate into the neighbouring cells, thus allowing 
energy transfer through collisions among neighbouring particles. The likelihood of 
such collision events can be controlled by the geometry of the cell. 



We assume that the scale separation, uq -C i/\y, is achieved. This involves specific 
conditions which depend on the geometry of the confining cells, but which will not 
concern us here. This assumption is to say that individual particles typically perform 
many collisions with the solid matrix of their confining cells, rattling about their cages 
at higher frequency than that of binary collisions. As will be shown in the next section, 
in this regime, Liouville's equation governing the time evolution of phase-space densities 
reduces to a master equation for the time evolution of local energies. The scale separation 
between the two collision frequencies v-q and z/w is the only parameter that controls the 
validity of this reduction. Moreover it becomes exact in the limit of vanishing binary 
collision frequency. 



4. Statistical evolution 



The phase-space probability density of the mechanical systems of N confined hard 
spheres is specified by the distribution function Pn(vi, Vi, . . . , rjy, where r a and 
v a , a = 1, . . . , N, denote the ath particle position and velocity vectors. The index a 
stands for the label of the confining cells. For our system, as is customary for hard 
sphere dynamics, this distribution satisfies a pseudo-Liouville equation [H], which is well 
defined in spite of the singularity of the hard-core interactions. This equation, which 
describes the time evolution of p^ is composed of three types of terms: (i) the advection 
terms, which account for the displacement of the moving particles within their respective 
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billiard cells; (ii) the wall collision terms, which account for the wall collision events, 
between the moving particles and the solid matrix of their confining cells; and (iii) 
the binary collision terms, which account for binary collision events, between moving 
particles belonging to neighbouring cells : 

N -I N 

dtPN = E [~«« ■ d ra + K^] p N + - E B^p N . (11) 

a=l ^ 0,6=1 

The first two terms account for the motion of individual particles within their respective 
cells, whether free advection or wall collisions, and the third term is the binary collision 
operator which, written in terms of the relative positions r ab and velocities v ab of 
particles a and b, and the unit vector e ab that connects them, is 

B {a > b) p N (. . . , r , « a , . . . , r b , v bl ...) = 

(2p m ) 2 / de ab (e ab ■ v ab ) \5(r ab - 2p m e ab ) 

xp N ( ...,r a ,v a - e ab (e ab ■ v ab ), ...,r b ,v b + e ab (e ab ■ v ab ), . . . ) 

-5(r ab + 2p m e ab )p N ( . . . , r a , v a , . . . , r b , v b , . . . )] . (12) 

4-1. Equilibrium states 



Equilibrium states of equation (11) are product measures such as a Maxwellian 



distribution with inverse temperature /?, 



:'..Y 
2 



(can)/ \ 1 / m f3 \ " -!2±3( v 2 + +v 2 ) , > 

p K N >(r 1 ,v 1 ,...,r N ,v N ) = I — 1 e = m+ '" + "» J , (13) 

where |£(AT)| denotes the volume accessible to all iV particles, or micro-canonical 
measures of totally energy E, 

(mic)/ v _ 1 l^-J-fc ( m \ 2 ^ ( m 



"••">=iwS Hi? 1 --^' (14) 



where the Gamma function is 

3iV-2 

3AT-2 



/3iV\ _ (3iV-2)!!V27r2-— , iVodd, 
1 I 2 J " (3iV-2)!!2-^, iV even. (15j 



4-2. Local equilibrium states 

Under the assumption that particles typically collide more often with the solid matrix 



of their confining cells than among nearest neighbours, the time-evolution (11) is 
dominated by the local terms so that the distribution p N rapidly evolves toward a local 
equilibrium distribution, which depends upon the energy variables {ei, . . . , e^} only. 

We denote by P/v(ei, . . . , ejv, t) this object, invariant under the local evolution terms, 
defined in terms of the phase-space distribution function of the mechanical systems of 
iV confined particles Pn(ti, vi, . . . , rjy, VN,t) according to 

P N (e 1} ...,e N ,t) = J f[ dr a dv a 5 (e a - p N (r 1 ,v 1 , . . . ,r N , v N ,t) . (16) 
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Notice that the normalization of pn implies that of Pn- 

In terms of the energy distributions, the corresponding equilibrium measures are 

3N 

y/ex ■ ■ ■ e N e'^ + - +€N) . (17) 




p(mic)/ x 



N 



One can check that both these measures are normalized. We notice that the one- 



and two-point distribution of the canonical and micro-canonical measure (18) have the 
forms (N > 3) 

Pt n) ^) = 2 -^V~ee-^, (19) 



Pr\t)= " .^a/' ' l-~ , (21) 



Pt n \^ e b ) = ^ y/e^-e^ , (20) 

7T 

,(mic), , _ _J_ r(|iV) / e x 

^lr[|(iv-i)] ve l £ 

. „ / 3 Ar x 3iV-8 

p (mic) r 4r(^iV) / e a e b \ 2 

(e., 6 6 ) - ^srdjy _ 3 ) ^ y 1 ~ E " £ J • (22) 

The canonical ensemble marginals are straightforward. As of the micro-canonical 
ensemble, the marginals are obtained from the iV-point distribution (18) by integrating 
out iV — 2 of the N energy variables. The integration over the first of N — 2 variables 
takes care of the delta function. Each subsequent integration has the form 

3 



dy^(x-y) a ~x' +a , (23) 

where x = 1 — e\/N — ... — e n /N. The successive values of a are 

\ when integrating over ejv-i > 

2 when integrating over e,/v-2 , 

\ when integrating over e7v-3, 

5 when integrating over ejv-4, 



and so on until we obtain (22) after N — 2 integrations. 

Letting E = 3/2N/3~~ 1 into (21 ), and taking the limit N — > oo, we recover 



,3 



lim P ( r\e) = 2 -^V-ee-^, (25) 



which is the one-particle distribution of the canonical measure (17). 
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The energy moments with respect to the one cell distributions ( 19 ), for the canonical 



ensemble, and (21), for the micro-canonical ensemble, are easily computed. Let (.} ( 



and (.) m ic denote the averages with respect to the distributions (19) and (21) at energy 
E = 3/2N(3~ 1 respectively. The corresponding nth moments of the energy are 



2 r( 3 +n 

7T \2 



(2n + l)H 

2™ ' 
J_ (3N\ n r(f +n)r(fAQ 

(2n + l)!!(f) n r(|iV) 



(26) 



2 n T(±N + n) 
In particular the first few moments are 

(/3e)mic = r , 



(27) 



15 3N 
~J3N + 2 ' 
105 9N 2 



(2f 



8 (3iV + 4)(3iV + 2) ' 

Thus only the first one is independent of N and equal to its canonical counterpart. The 
remaining ones are only equal to their canonical counterparts up to 0(1 /N) corrections. 



4-3. Local equilibrium closure 



Taking the time-derivative of equation (16) and substituting equation (11), only the 



term involving energy exchanges between neighboring particles contribute : 
d t P N {e u . . . ,e N ,t) = J JJdrid«i«J (e { - ^vfj d t p N (r u v 1} . . . , r N , v N , t) , 



1 N N , 

2 J II drid« f (J U - 

a, 6=1 i=l 



~v*) B^p N (r 1 ,v 1 ,...,r N ,v N ,t){29) 



The action of B^ a ' b ^ on p^ is given according to equation (12). In order to obtain a 



closed equation for P N , we assume that the distribution p N is a locally micro-canonical 
one, i.e. it depends only on the energy distribution and is thus given according to 

,N , N 



pjv(n,«i, • • -,r N ,v N} t) 



m 



lir) N \£(N)\ 



i=l 




P iV (ei,...,e^,t).(30) 



The other prefactors are so chosen that this measure is normalized, viz. 

N 

I JJ drjdv^iv(r*i, v x , . . . , r N , v N , t) = 1 . 

i=i 



(31) 
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5. Master equation 



The validity of equation (30) depends on the scale separation between the wall and 



binary collision frequencies which will be assumed throughout this article. This allows 
us to focus on the stochastic reduction of the energy exchange dynamics, considering the 
mesoscopic level description of the time evolution of probability densities as a stochastic 
evolution. 



Indeed, substituting equation (30) into (12), we write the time evolution of local 



equilibrium distributions (29) in the form of the master equation (JTl), whose kernel W 
can be identified from the computations above. The time evolution is thus specified by 
a master equation which accounts for the energy exchanges between neighbouring cells 
and makes no further reference to the collisional dynamics of confined hard spheres. 



5.1. Derivation of the master equation 



For each pair (a, b) of neighbouring cells, we have a contribution to equation (29) of the 
form 

jf[dndVi6 (e t - jv^ B^% N ({r t ,v t },t) 

N 

= 4 Pm / II dr i dv i / de afe (e ab • v ab )S ( e { - — t 

J i= i Je ab -v ab >0 V Z 



(32) 



x [5(r ab - 2p m e ab )p N ( . . . , v' a , v' b , . . . ) - 5(r ab + 2p m e ab )p N ( . . . , v a , v b , . . . )] 



where v' a = v a - e ab (e ab ■ v ab ) and v' b = v b + e ab (e ab ■ v ab ). Substituting equation ( |30D 
into this expression, we have the contributions 

, jV „ N 



ip 



m 



\C(N)\ 



/ YldridViS Ui 

J 8=1 ^ 



^rvf] j de ab (e ab ■ v ab ) 



(33) 



i=i 




x 5(r ab - 2p m e ab ) J f[ de^ x 8 | v • - | P N {ei, ...,e N ,t) 

N 

-5(r ab + 2p m e ab ) J J[ de i e l 1 5 
where v[ = Vi for i ^ a, b and 



m 

r, - M— I fjv(ei, . . .,e N ,t) 
m 




v 2 a - {e ab ■ v a ) 2 + (e ab ■ v b ) 2 , 



(34) 



j b - \l v t + K e ab ■ v a y - {e ab ■ v b y . 

Inserting factors J dq5{rj ± m/2[(e ab ■ v a ) 2 — (e ab ■ v b ) 2 }) in both of the terms between 



brackets, we write e[ — i ^ a, b, and e' a — e a + 77, e' b = e b — rj. The contribution (33) 
thus transforms into 



m 



(87r) 7V |£(A^)| V2 



m\ 2 ' 



i=i 



m , 



e a b-v ab >0 



de ab (e ab ■ v ab ) 
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x / jjde i e i 1 (y( > /e' i - > /e i )Pj V (ei,...,ejv,t) 
J i=i 

dr}8 [v - — \(e ab ■ v a ) 2 - (e ab ■ v b 

N 
i=l 

(8ny v \L,(N)\ J V 2 / Je ab -v ab >o 

x S(r ab - 2p m e ab ) J di]5 [n + — [(e ab ■ v a ) 2 - (e ab ■ v b 



ab( e ab ■ V ab ) 



X 



e\ . . . e 



=-P/v(e'i, • • • , %, t) 



N 



dr]5 [n - — \{e ab ■ v a ) 2 - (e ab ■ v b ) 



x 



-Pn(^i, ■ ■ ■ , ejV) t) 



. - , , (35) 

V e l • • • e N 

where we used the identity de/e = 2dy / e/y / e and carried out the e integrations. We 
proceed by performing all the and Vi integrations but for i — a,b. For the latter, we 
change variables from (r a ,r b ) to (R ab ,r ab ), the center of mass and relative coordinates 
respectively. The latter integration can be carried out, with outcome r ab = ±2p m e ab . 
We are thus led to contributions 

2 S 

8tt 2 |£(2)| 



m 



dR ab dv a dv b 5 e a - —v a )8[e b 



m 



de ab (e ab ■ v ab ) 



x 



1 



drj5 f T) + — [(e ab ■ v a ) 2 - (e ab ■ v b ) 



xP N (e 1 ,...,e a + r],e b -r],...,e N ,t) 

^= f dr]5 (rj - ^-[{e ab ■ v a ) 2 - (e ab ■ v b ) 2 ]) P N (ei, ...,e N ,t) . (36) 

yje a e b J \ l J 

Notice that we have made the assumption that the position integrals of the iV — 2 
particles not involved in the collision factorise^J 

\C(N)\ = \C(2)\\C(N - 2)\ . (37) 

Though this is an approximation, it becomes exact under the assumption of scale 
separation, which implies |£(iV)| = 1/2(1)1^. Comparing the above expression to 
equation Q, we identify the kernels associated with the gain and loss terms, 



W{e a + r],e b -r]\e a ,e b ) 



Pm m 



dQdR 



8tt 2 |£(2)| 

All the cells are taken to be identical throughout. 
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= / dv a dv b (e ab ■ v ab )5 (e a 

{e a + r])(e b -ri) ^ a t-v ab >o \ 

dttdR 



m 



8^\C(2)\ 
1 



x 



e a + r}){e b -r}) J ^ b -v ab <o 



dv' a dv' b (e ab ■ v' ab )5 ( e a + 77 - 2 



x5 ( e 6 — 77 — ^-v 1 ^) 8(r]+ '^r[{e ab -v 



m , 



.' ^2 



W(e a ,e b \e a -r),e b + 77) 



2 

Pm m 

ln 2 \C(2)\ 



2 

dVldR 



{Sab ■ V' a ) 2 ) 



(3* 



x 



/ dv a di; 6 (e a;) • « a6 )<5 (e a 

Je ah -v nh >0 V 



III 



\ftatb Je ab -v ab >0 



2 a 



ni 



m , 



x5[e b - —v b )5 [r] - —[{e ab ■ v a ) - (e ab ■ v b 



(39) 



where, in the expression of the gain term (38), we changed the post-collisional velocity 
variables (v a ,v b ) to the pre-collisional velocities (v' a ,v' b ). 

Without loss of generality, we assume v a and v b are measured with respect to 
referential whose 2-axis point along the unit vector e ab . In such a case, we have 
e a b • v a = v a cos 4> a and e ab ■ v b — v b cos <p b . Thus the volume integral / dfidi? decouples 
from the velocity integrals which yield 



e ab -v' ab <0 



dv' a dv' b (e ab ■ v' ab )5 ( t a + 77 - ) 5 ( e b - 77 



m ,2 



(40) 



m , 



xS [V + y[(e a b ■ v b ) -(e a b-v' a 



I \2l 



71 



2\l 



in 



~a c b 



dx a dx b (J e' b x b - Je' a x a )S(r] + [e b x b - e a x a \) , 



/ dv a dv b (e ab ■ v ab )5 (e a - ^-v 2 a ) 5 (e b 

Je ab -v ab >0 \ A J \ 



m 



(41) 



x8 ( 77 - y[(e ab • v a ) 2 - (e ab • -y 6 ) 2 ; 



7T 



2\i 



r?) 



/ 



\/^iX a >^/I b ~X b 



dx a dx b (^fe^x a - y^x^Sirj - [e a x 2 a - e b x 2 b }) 



Notice that the two integrals are identical upon exchanging the roles of a and h and 
77 — > —77, except for the energies e' a = e a + 77 and e' b = e b — 77 in the expression (40). Thus 
the kernels (38)-(39) have expressions 

W(e a + r),e b -r)\e a ,e b ) ~ 1 



1 2 
in 



dttdR 



x 



e'x b >We' a x 



m\C(2)\ 

dx a dx b (J7 b x b - J7 a x a )5(r] + [e' b x 2 b - e' a x 2 a }) 



(42) 
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W(e a ,e b \e a -r],e b + rj) 



i 2 
in 



dQdR 



m\C(2)\ 

dx a dx b (y/e^x a - yJF h x b )8(r) - [e a x 2 a - e b xl\) 



x 

Notice the symmetry between the kernels of the gain and loss terms, i. e. 

W{e a + r],e b - r]\e a , e b ) = W{e' a , e' b \e' a - rj, e' b + rj) . 
This is however no to say that the kernels are symmetric. In fact 

W(e a + r],e b - rj\e a , e b ) ^ W(e a , e b \e a + rj, e b - rj) , 



(43) 

(44) 
(45) 



2 p 



dQdR< 



which can be seen by rewriting equation (|42|) directly in terms of e a and e b 
W(e a + r],e b -r)\e a ,e b ) 

L 



x 



The symmetry is therefore the following 



m\£(2)\J - "V (e a + rj)(e b - V ) 
dx a dx b (^fe^x a - ^fe b x b )5(ri + [e a x 2 a - e b xl\) . 



(46) 



(e Q + rj)(e b - r])W(e a +r], e b -r]\e a} e b ) = ^J7^~ b W{e a , e b \e a +r], e b -rj) , (47) 
consistent with the form of the equilibrium distributions (Jl7|). This is to say that, 



despite the lack of symmetry (45), detailed balance is recovered: 



W{e a +r),e b - r)\e a , e fe )P (cq) (. . . , e a + rj, e b - rj, . . .) 
= W(e a , e b \e a + V ,e b - 7?)P (eq) (. ••,£«, e b , . . .) . 
Therefore d t P^ eq '(ei, . . . ejy) = 0, as must be. 

5.2. Computation of the loss term 
Going back to equation ( [43] ) , we can write 
8(V ~ e aX 2 a + e b x 2 b ) 



(4* 



(49) 



1 



2y/e a (ri + e b xl) 



'rj + e b x£ 



+ S\ x a + 



I rj + e b xi 



The validity of this expression requires rj + e b x b > 0, which is satisfied whenever rj > 
or 



abs(xb) > J — ^ , if rj < . 



(50) 



Provided this condition is fulfilled, we can carry out the x a -integration in equation (43). 



Considering the arguments of the delta functions in equation (49), the result of the 
integration would be trivial unless 



rj + t b x\ < e a & abs(x b ) < 



V 



e b 



(51) 
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The bounds of the avintegral are respectively 



x m i n = max[-l,- v /(e a -^)/e 6 ], 

(52) 



x max = min[l, yj (e„ - rf)/e b ]. 
which is 

x __ x . _ j l ) v«*-*» (53) 

■^max A min I e a -r) ^ _ 

£!) > '/ fc « fc & ' 

77 > 0: The condition on the integration domain transposes into the Heaviside 
step function with arguments ±1/77 + e b x b — y^x?,. Only the positive sign contributes 



to the kernel (43). Leaving aside the prefactors, we have 

1 f Xmax j fi y^bXb 



| C ^ max / 

W(e a ,e b \e a — 77, e 6 + 77) ~ — — / dar 6 1 



1 



where we used the symmetry x min = — x max . Using equation (53), we have 



f 7=, 0<77<e a -e fe (ife a >e 6 ) 



(55) 



5.2.2. 77 < 0: In this case, the ^-integral in (43) splits into two separate integrals, the 



first from x min , as defined above, to — x mic i, and the second from £ mid to x max , with 



x mid = J — - . (56) 
V e& 

Note that the bounds are well ordered : we always have x mi d < x max . 

The terms contributing to the kernel are ±\Jrj + e b x b — ^[t b x b for x b < 0. We have 

u/f 1 , ^ 1 f' Xinid a ( V^bXb 

W(e a ,e b \e a -rj,e b + r]) ~ — — / dx b \-2—f 

(57) 



V + £bX bl 

j V + e ^max I V + £ Ad 

V e a e b V e a e& 

Substituting the expressions of x max and x m id, we find 

f -7= , e a - e b < 77 < (if e 6 > e a ) , 
W(e a ,e 6 |e a - 77, 6^ + 77)- ^ ' n ( 58 ) 

5.2.3. To summarize: We can write for the loss term: 

W(e a , e b \e a - 77, e b + 77) = \[^T^ J <^dR (59) 
' max(0,e a - o,) < // < <„ . 



1 



\J max(e„,e t ) 



min(0, e a - e b ) < 77 < max(0, e a - e 6 ) 



-e& < 77 < min(0, e a - e b ) . 
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It is readily checked that W is symmetric under exchanging e a <-> e b and 77 — > —77, as it 
must be, and has the inverse units of time times energy. 

5.3. Gain term 



Going through the same steps as in Section 5.2 above, we find for the gain term: 



W(e a + rj,e b -r]\e a ,e b ) = / dftdi* (GO) 



m\C(2)\ 



1 



max(e b - e a , 0) < 77 < e 6 



(^feC-T?) ' min ( e & - e a , 0) < 77 < max(e 6 - e B , 0) , 
1 -e a < 77 < min(e fc - e a , 0) . 



Comparing with equation (59), we check that the loss and gain terms have the 



symmetry of equation (48). 



6. Collision frequency and heat conductivity 

We can rescale time by a factor (4p m ) 2 /(^/m^T\C(2)\) J dQdR, which amounts to 
converting the time units to units of the inverse square root of energy. The expression 
of the kernel thus simplifies to: 

W(e a ,e b \e a - r],e b + r]) = (61) 

_ e6 < v < min(0, e a - e b ) , 



/^xJ 7 = min(0, e a - e b ) < 77 < max(0, e a - e b ) , 

V 8 I ymax(e a ,e 6 ) 

max(0, e a - e b ) < 77 < e a . 



A graphical representation of this function is displayed in figure |2j 
The rate of energy exchange ^ is 

27 [ e b >e a 



"M = V X - ' (62) 



with equilibrium average taken with respect to the measure (17), 

u B = -^g J de a de b u(e a , e b )^T b e-^ +eb)/T = Vf . (63) 
The micro-canonical average of the collision frequency is the quantity 

u B (N) = J de a de b u(e a , e b )P^ ic \e a , e b ) = Vf[l + 0(1/N)} . (64) 

The amount of heat transfer is 

27. . ( e b >e a 



i(^e 6 ) = ~(e tt -e 6 )x<{ 3 V ?£6 ' "/ (65) 



which is equivalent to 



j( e a, e b ) = -(e a - e b )z/(e a , e 6 ). (66) 




Thus, with the help of the identity 

i J de a de b u(e a ,e b )(e a -e b ) 2 ^T b e-^ + ^ T = Vf, (67) 

we can compute the average of the heat current weighted by (e a — e b )/2 with respect to 
local thermal equilibrium measure to obtain the heat conductivity: 
2 



7T 



y5 



J de a de b j(e a , e b )(e a - t b )^/e^e b e 



(e a +e b )/T 
(e a +e b )/T 



J de a de b u(e a ,e b )(e a - e b ) 2 y/e^e~ b e~ 
/ f, (6S 



The ratio between the heat conductivity and collision frequency is therefore unity, as 
announced 

The same result is obtained from the equilibrium average of the second moment of 
the heat transfer (JVJ) : 

fK~ ( 35el(e b -e a )+2U b e 2 a +lle 3 a 

hi \ — ^ J ' b — a ' ^«n\ 

^ 6a ' ^ ~ 420 I ^l(e a -e b )+2U a el+Uel l bJ ) 
I ' a ~ 6 ' 

which yields an alternative derivation of the heat conductivity by taking the equilibrium 
average of this quantity: 
2 



TIT 5 



y de a de 6 /?,(e a ,e fe ) v /e a e 6 e" 



(e a +e b )/T 



Vf. (70) 
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This result is universal for three-dimensional systems in the sense that it holds for 
confined particles interacting through hardcore collisions, independent of the shape of 
the confining cells and temperature, at least so long as the local dynamics is ergodic on 
the constant energy surface. 



7. Numerical scheme 



Equation (J5]) can be verified by direct numerical computations of the master equation 
PJ) according to the scheme described in [B], based on Gillespie's algorithm [HJ HO] . 
We recall that the Monte-Carlo step necessitates two random trials. The first random 
number determines the time that will elapse until the next energy exchange event. The 
second one determines which one out of all the possible pairs of cells will perform an 
exchange of energy and how much energy will be exchanged between them. 

Thus let £ and \ be these two random numbers, uniformly distributed on the unit 
interval and let us consider an array -whether one-, two- or three-dimensional- of N 



pain 



of energy variables {e^, e^j. . 
The first random number, £, determines the time to the next energy exchange event, 
denoted r, according to the Poisson distribution whose relaxation rate is given by the 
sum of the respective exchange rates u(e[ , e% j between pairs of energies {e^ , ejj }, 

The second random number, x, determines the pair n of energies (ei , j undergoing 
the energy exchange and how much energy rj they exchange, according to 

-y "n— 1 

X 



EM*® 4°) 



E "(4°, 4°) + f\ n) */W(e<?\ 4 n) |cS B) - V, 4 n) + i? 



which must be solved for 77. Thus, letting 



x 



i=l 



i=0 



and assuming 



n) (n) 
) e 2 



we have 
3 



V 



< x - d V 'w(e ( r\ e^\et ] ~ V, ^ +v)< 



(72) 

(73) 
(74) 



2^4"'4" )\ 1/3 



>) 
~ e 2 ) 



o Pi ( / W / 2 • I (n) 

ZxJ^maxi ye\ , y e 2 J — |mm( , e 

(n) 



(n) » 
2 



>) 
e 2 



+min(0, €1 

o f (") f («)\ 1 / 3 



n) r (n) 
2 



>) 
e l 



3vr 



2/3 



^(ef^ej ) <x<v 2 (e^,e 

v 2 \£\ , £2 



(75) 



, (n) (n) 

< x < v( e\ , e\' 



|| The corresponding number of energy variables depends on the dimensionality of the lattice. In 
dimension 1, each energy variable is involved in two pairs; in dimension 2, four pairs; in dimension 3, 
6 pairs. 
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In the above equation, we introduced the intermediate bounds v\ and v% which are 
computed through the partial integrals 

f d V 'W(e a , e b \e a - rf, t b + rf) = (76) 

J — e h 



X < 



2min(e a ,ei,)+3[T?-min(0,£ a -£j,)] 
6max(^i^ )V ^) 



e a +e t +2max(e a ,e t ) _ (e a -r))^ 



-e b < 7] < min(0, e a - e b ) , 
min(0, e a - e b ) < rj < 

max(0,e a - e h ) 

max(0, e a - e b ) < r] < e a . 



In particular, 



vi(£a,e b 



min(0,€ a -et) 



dr)W(e a ,e b \e a -f],e b + r]) 



2tt min(e a , e b ) 
6 max( v /e^, ^/F b ) ' 

max(0,e a -e 6 ) 

dr]W(e a ,e b \e a - r],e b + i]) 



(77) 



£6 



y/^r min(e a ,e fe ) 
6 max( > /i^, ^/e~ b ) ' 
= i/(e OJ e 6 ) - vi{e a ,e b ) . (78) 

We remark that the scheme is explicit here and does not require using root finding 
algorithms to determine the amount r] of energy exchanged, as was the case with the 
kernel associated with two-dimensional underlying dynamics in |6]. This considerably 
reduces the amount of CPU time necessary to implement the algorithm. 



7.1. Equilibrium simulation 

We briefly discuss the results of equilibrium simulations performed along this scheme 
on one- and two-dimensional arrays of iV cells^Jwith periodic boundary conditions, i.e. 
identifying cells iV + 1 and 1 in the one-dimensional case and working with a square 
lattice of size y/N x \J~N and identifying the first and last columns and rows in the 
two-dimensional case. 

The top panel of figure [3] shows the results of a numerical computation of the first 
three energy moments for the one-dimensional lattice with different values of iV up to 
iV = 250 and compares it to equation (28 ). The accuracy of the calculation demonstrates 



the invariance of the one-cell distributions (21), and by extension, of the micro-canonical 



distribution (18). This is further validated by a computation of the energy exchange 



frequency, which involves the two-cell distribution (22). The bottom panel of figure 3 



thus shows the results of a numerical computation of this quantity for values of N up 
to iV = 250, together with a comparison to a numerical evaluation of equation (64). 
Similar results are obtained for two-dimensional lattices. 

% Notice the change of notation: N refers to the system size here and not the number of pairs of 
neighbouring cells. 
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Figure 3. (Top) Energy moments, ((e/T)™), for n — 1,2,3, for values of N up 
to N — 250. The blue dots are the results of numerical computations of these 
quantities, using the scheme described above, and the solid lines correspond to their 



theoretical values derived in equation ( 28 1 . The dashed line correspond to the canonical 
expectations (26 1. (Bottom) Micro-canonical energy exchange frequency vs. N. The 



blue dots are the results of numerical computations of the energy exchange frequency 
for the corresponding value of N. The solid red line is the numerical integration of the 



equation (64) 



The computation of the mean-squared displacement of the Helfand moment, 
H(t) = Y,a=i ae a(t), yields the conductivity in the limit of large system sizes, 

k = lim - \ NO lim / — AH(T n ) 2 ) . (79) 
As seen in figure |4l the conductivity is identical to the collision frequency, in accordance 
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with equation ([5]): 

k ( 0.9994 ±0.0016, (ID), 
ub ~ { 0.9993 ±0.0026, (2D). 
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Figure 4. Ratio between the mean squared displacement of the Helfand moment 
and energy exchange frequency for one- (green) and two-dimensional (blue) lattices 
of cells vs. the size of the lattices N. The ratio k/vb is the infinite N extrapolation 
of these data, which we evaluate by linear regression (dashed red line). Fitting the 
data points corresponding to N > 9 (of which only a fraction are displayed here) 
with weights inversely proportional to the sizes of their error bars, we obtain, for 
the one-dimensional lattice, k/v-q = 0.9994 ± 0.0016, with N up to 400, and, for the 
two-dimensional lattice, k/v b = 0.9993 ± 0.0026, with N up to 400 (= 20 x 20). 



7.2. Non- equilibrium simulation 

Non-equilibrium boundary conditions are easily implemented on one-dimensional 
lattices by thermalising the boundary cells at respective temperatures T_ = 0.5 and 
T + = 1.5 according to the scheme detailed in [6j. 

The heat conductivity is evaluated by computing the stationary heat current, Jh, 
and comparing it to the temperature gradient according to 

T — T 

Jk = -k^ w ^, (81) 
in the limit of large system sizes N. 

(2) 

Denoting by P a a+1 (e,e') the two-point marginal of the stationary state associated 
with cells a and a + 1 at respective energies e and e', the heat current is 

j n = j dede'j(e,e')P! l % 1 (e,e'). (82) 

A numerical computation of this quantity was performed for different system sizes, 
as shown on the left panel of figure [5| The result of the infinite N extrapolation yields 

— = 1.00016 ±0.0005. (83) 
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According to Fourier's law, the corresponding local temperature profile is expected 



to be 



T 



1 3 3 fi 3 



which is confirmed numerically as seen in the right panel of figure [5] 
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n/(N+l) 
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Figure 5. (Left) Ratio between the heat current (82) and local temperature gradient 
under thermal boundary conditions at the temperatures T + = 1.5 and T_ = 0.5, 
corresponding to the temperature profiles shown above. The system sizes N here 
range from N — 1 to N = 100. The ratio k/v-q is computed from the infinite N 
extrapolation of these data, evaluated by linear regression with respect to 1 /N (solid 
red line). Fitting the data points weighed according to the sizes of their error bars, we 
obtain k/v-q — 1.00016 ± 510~ 4 . (Right) Temperature profiles obtained in the non- 
equilibrium stationary states of systems in contact with stochastically thermalised cells 
at n = (T_ = 1/2) and n = N + 1 (7+ = 3/2). The thick black line corresponds 
to the profile ( p4| expected from Fourier's law. Different system sizes are displayed, 
going from N = 5 to N = 100. The data is barely distinguishable from the curve (84 1 
when N is sufficiently large. 



8. Conclusions and perspectives 

To summarize, we have successfully extended the results presented in [B] to the confining 
dynamics of rarely interacting hard spheres. Our results thus further validate the claim 
that the identity between heat conductivity and collision frequency in such systems of 
confined particles interacting through rare hard core collisions is largely independent of 
the details of the underlying dynamics. 

In a forthcoming publication, we will show that the dimensionality of the underlying 
dynamics is indeed arbitrary. In particular, one can consider systems which consist of a 
solid structure of confining pores of arbitrary shapes in which gas particles of arbitrary 
numbers are trapped. The geometry of the pores must be such that each isolated pore 
contains a fixed number of hard spheres with micro-canonical equilibrium measure whose 
energy is the sum of the kinetic energies of the gas particles. Under the assumption that 
relaxation to this local equilibrium measure takes place on time scales smaller than the 
time scales of energy exchanges between neighbouring pores, a master equation of the 
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form ([T]) can be derived to describe the energy exchange process. The corresponding 
stochastic kernel will in this case depend on the geometry of the pores involved in the 
energy exchange and, in particular, on the precise number of particles in those pores. 
The symmetry relations (|8j) will thus reflect local properties of the system whose overall 
average yields the corresponding macroscopic properties. In this sense, one expects that 
Fourier's law can be derived in a rather general setting for such systems of confined 
particles, and the value of the heat conductivity expressed in terms of the frequency of 
energy exchanges between the cells. 
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